二維非穩(wěn)態(tài)導(dǎo)熱隱式格式求解Matlab代碼
%%感謝up主 @計(jì)算傳熱學(xué)大叔 無(wú)敵的教學(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=6000;
x=3;
y=2;
dt=10;
dx=0.1;
dy=0.1;
sizet=t/dt;
sizex=x/dx;
sizey=y/dy;
%此處不能像1D那樣投機(jī)取巧——廣義源項(xiàng)b必須是一個(gè)列向量才能求解方程，原有的三維數(shù)組/暴力求解方法失效
%可以試著降維成2維矩陣——對(duì)每個(gè)時(shí)間步，每一"行（x）"網(wǎng)格對(duì)應(yīng)一個(gè)對(duì)角陣，對(duì)角矩陣沿對(duì)角線連接。
%這樣做的好處是可以將整個(gè)計(jì)算域內(nèi)的所有網(wǎng)格溫度一起算，特別適合迭代求解。
Tmatrix=zeros(sizet,sizey,sizex);%時(shí)間，行，列
for i=1:1
? ?Tmatrix(i,:,:)=T0;
end
%% 計(jì)算系數(shù)矩陣和源項(xiàng)矩陣
AP1=zeros(sizey,sizex);%ap1系數(shù)矩陣
AE1=zeros(sizey,sizex);%ae1系數(shù)矩陣
AW1=zeros(sizey,sizex);%aw1系數(shù)矩陣
AN1=zeros(sizey,sizex);%an1系數(shù)矩陣
AS1=zeros(sizey,sizex);%as1系數(shù)矩陣
%沒(méi)設(shè)ap0，直接在b向量里計(jì)算就行
B=zeros(sizey,sizex);%廣義源項(xiàng)矩陣
for j=1:1%按列去分
? ?for r=1:1 %矩形區(qū)域左下頂點(diǎn)，dxw=0.5dx,dys=0.5dy
? ? ? ?ap0=dy*dx*rho*cp/dt;
? ? ? ?AE1(r,j)=k*dy/dx;
? ? ? ?AW1(r,j)=2*k*dy/dx;
? ? ? ?AS1(r,j)=2*k*dx/dy;
? ? ? ?AN1(r,j)=k*dx/dy;
? ? ? ?AP1(r,j)=dx*dy*rho*cp/dt+k*dy/dx+2*k*dy/dx+k*dx/dy+2*k*dx/dy;
? ? ? ?%ap1 ? ?=dy*dx*rho*cp/dt+k*dy/dxe+k*dy/dxw+k*dx/dyn+k*dx/dys;
? ? ? ?%廣義源項(xiàng)需要用到上一時(shí)刻的溫度，所以放到時(shí)間步內(nèi)去算
? ?end
? ?for r=2:sizey-1 %矩形區(qū)域左邊，dxw=0.5dx
? ? ? ?%不用算ap0了，因?yàn)槎家粯?/p>
? ? ? ?AE1(r,j)=k*dy/dx;
? ? ? ?AW1(r,j)=2*k*dy/dx;
? ? ? ?AS1(r,j)=k*dx/dy;
? ? ? ?AN1(r,j)=k*dx/dy;
? ? ? ?AP1(r,j)=dx*dy*rho*cp/dt+k*dy/dx+2*k*dy/dx+k*dx/dy+k*dx/dy;
? ? ? ?%ap1 ? ?=dy*dx*rho*cp/dt+k*dy/dxe+k*dy/dxw+k*dx/dyn+k*dx/dys;
? ?end
? ?for r=sizey:sizey %矩形區(qū)域左上頂點(diǎn)，dxw=0.5dx,dyn=0.5dy
? ? ? ?%不用算ap0了，因?yàn)槎家粯?/p>
? ? ? ?AE1(r,j)=k*dy/dx;
? ? ? ?AW1(r,j)=2*k*dy/dx;
? ? ? ?AS1(r,j)=k*dx/dy;
? ? ? ?AN1(r,j)=2*k*dx/dy;
? ? ? ?AP1(r,j)=dx*dy*rho*cp/dt+k*dy/dx+2*k*dy/dx+2*k*dx/dy+k*dx/dy;
? ? ? ?%ap1 ? ?=dy*dx*rho*cp/dt+k*dy/dxe+k*dy/dxw+k*dx/dyn+k*dx/dys;
? ?end
end
for j=2:sizex-1
? ?for r=1:1 %矩形區(qū)域中間下邊，dys=0.5dy
? ? ? ?AE1(r,j)=k*dy/dx;
? ? ? ?AW1(r,j)=k*dy/dx;
? ? ? ?AS1(r,j)=2*k*dx/dy;
? ? ? ?AN1(r,j)=k*dx/dy;
? ? ? ?AP1(r,j)=dx*dy*rho*cp/dt+k*dy/dx+k*dy/dx+k*dx/dy+2*k*dx/dy;
? ? ? ?%ap1 ? ?=dy*dx*rho*cp/dt+k*dy/dxe+k*dy/dxw+k*dx/dyn+k*dx/dys;
? ?end
? ?for r=2:sizey-1 %矩形中間區(qū)域中間，主體部分
? ? ? ?AE1(r,j)=k*dy/dx;
? ? ? ?AW1(r,j)=k*dy/dx;
? ? ? ?AS1(r,j)=k*dx/dy;
? ? ? ?AN1(r,j)=k*dx/dy;
? ? ? ?AP1(r,j)=dx*dy*rho*cp/dt+k*dy/dx+k*dy/dx+k*dx/dy+k*dx/dy;
? ? ? ?%ap1 ? ?=dy*dx*rho*cp/dt+k*dy/dxe+k*dy/dxw+k*dx/dyn+k*dx/dys;
? ?end
? ?for r=sizey:sizey %矩形區(qū)域中間上邊，dyn=0.5dy
? ? ? ?AE1(r,j)=k*dy/dx;
? ? ? ?AW1(r,j)=k*dy/dx;
? ? ? ?AS1(r,j)=k*dx/dy;
? ? ? ?AN1(r,j)=2*k*dx/dy;
? ? ? ?AP1(r,j)=dx*dy*rho*cp/dt+k*dy/dx+k*dy/dx+2*k*dx/dy+k*dx/dy;
? ? ? ?%ap1 ? ?=dy*dx*rho*cp/dt+k*dy/dxe+k*dy/dxw+k*dx/dyn+k*dx/dys;
? ?end
end
for j=sizex:sizex
? ?for r=1:1 %矩形區(qū)域右下頂點(diǎn)，dxe=0.5dx,dys=0.5dy
? ? ? ?ap0=dy*dx*rho*cp/dt;
? ? ? ?AE1(r,j)=2*k*dy/dx;
? ? ? ?AW1(r,j)=k*dy/dx;
? ? ? ?AS1(r,j)=2*k*dx/dy;
? ? ? ?AN1(r,j)=k*dx/dy;
? ? ? ?AP1(r,j)=dx*dy*rho*cp/dt+2*k*dy/dx+k*dy/dx+k*dx/dy+2*k*dx/dy;
? ? ? ?%ap1 ? ?=dy*dx*rho*cp/dt+k*dy/dxe+k*dy/dxw+k*dx/dyn+k*dx/dys;
? ?end
? ?for r=2:sizey-1 %矩形區(qū)域右邊，dxe=0.5dx
? ? ? ?%不用算ap0了，因?yàn)槎家粯?/p>
? ? ? ?AE1(r,j)=2*k*dy/dx;
? ? ? ?AW1(r,j)=k*dy/dx;
? ? ? ?AS1(r,j)=k*dx/dy;
? ? ? ?AN1(r,j)=k*dx/dy;
? ? ? ?AP1(r,j)=dx*dy*rho*cp/dt+2*k*dy/dx+k*dy/dx+k*dx/dy+k*dx/dy;
? ? ? ?%ap1 ? ?=dy*dx*rho*cp/dt+k*dy/dxe+k*dy/dxw+k*dx/dyn+k*dx/dys;
? ?end
? ?for r=sizey:sizey %矩形區(qū)域右上頂點(diǎn)，dxw=0.5dx,dyn=0.5dy ? ? ? ?
? ? ? ?%不用算ap0了，因?yàn)槎家粯?/p>
? ? ? ?AE1(r,j)=2*k*dy/dx;
? ? ? ?AW1(r,j)=k*dy/dx;
? ? ? ?AS1(r,j)=k*dx/dy;
? ? ? ?AN1(r,j)=2*k*dx/dy;
? ? ? ?AP1(r,j)=dx*dy*rho*cp/dt+k*dy/dx+2*k*dy/dx+2*k*dx/dy+k*dx/dy;
? ? ? ?%ap1 ? ?=dy*dx*rho*cp/dt+k*dy/dxe+k*dy/dxw+k*dx/dyn+k*dx/dys;
? ?end
end
%% 開(kāi)始迭代
for i=2:sizet
? ?T1=T0+zeros(sizey,sizex);%迭代初值
? ?T2=T1;
? ?ErrTol=1e-6;
? ?Err=1;
? ?while Err>ErrTol
? ? ? ?for j=1:1%左側(cè)，整個(gè)這部分的源項(xiàng)都收納了TW*AW1一項(xiàng)，于是迭代過(guò)程不用寫(xiě)
? ? ? ? ? ?for r=1:1%左下角，南北是r，是行，從南到北是1到sizey；東西是j，是列，從西到東是1到sizex。
? ? ? ? ? ? ? ?B(r,j)=s*dx*dy+ap0*Tmatrix(i-1,r,j)+TW*AW1(r,j)+TS*AS1(r,j);
? ? ? ? ? ? ? ?%源項(xiàng)不隨溫度迭代而改變，這是個(gè)確定的數(shù)。
? ? ? ? ? ? ? ?T2(r,j)=(B(r,j)+AN1(r,j).*T1(r+1,j)+AE1(r,j)*T1(r,j+1))./AP1(r,j);
? ? ? ? ? ? ? ?%T2是迭代計(jì)算所得的溫度，Tmatrix是最終記錄到結(jié)果矩陣中的溫度
? ? ? ? ? ?end
? ? ? ? ? ?for r=2:sizey-1%左側(cè)中部
? ? ? ? ? ? ? ?B(r,j)=s*dx*dy+ap0*Tmatrix(i-1,r,j)+TW*AW1(r,j);
? ? ? ? ? ? ? ?T2(r,j)=(B(r,j)+AN1(r,j)*T1(r+1,j)+AS1(r,j)*T2(r-1,j)+AE1(r,j)*T1(r,j+1))./AP1(r,j);
? ? ? ? ? ? ? ?%源項(xiàng)比左下角少了一個(gè)TS*AS1，迭代項(xiàng)多了一個(gè)T2(j,r-1)*AS1，注意是T2，因?yàn)門(mén)2(j,r-1)已經(jīng)算出來(lái)了
? ? ? ? ? ?end
? ? ? ? ? ?for r=sizey:sizey
? ? ? ? ? ? ? ?B(r,j)=s*dx*dy+ap0*Tmatrix(i-1,r,j)+TW*AW1(r,j)+TN*AN1(r,j);
? ? ? ? ? ? ? ?T2(r,j)=(B(r,j)+AS1(r,j)*T2(r-1,j)+AE1(r,j)*T1(r,j+1))./AP1(r,j);
? ? ? ? ? ? ? ?%相比左側(cè)中部就是源項(xiàng)多了一個(gè)AN1*TN，然后迭代項(xiàng)少了。
? ? ? ? ? ?end
? ? ? ?end
? ? ? ?for j=2:sizex-1%中部
? ? ? ? ? ?for r=1:1%中部底邊
? ? ? ? ? ? ? ?B(r,j)=s*dx*dy+ap0*Tmatrix(i-1,r,j)+TS*AS1(r,j);
? ? ? ? ? ? ? ?T2(r,j)=(B(r,j)+AN1(r,j)*T1(r+1,j)+AE1(r,j)*T1(r,j+1)+AW1(r,j)*T2(r,j-1))./AP1(r,j);
? ? ? ? ? ?end
? ? ? ? ? ?for r=2:sizey-1%主體部分
? ? ? ? ? ? ? ?B(r,j)=s*dx*dy+ap0*Tmatrix(i-1,r,j);
? ? ? ? ? ? ? ?T2(r,j)=(B(r,j)+AN1(r,j)*T1(r+1,j)+AE1(r,j)*T1(r,j+1)+AW1(r,j)*T2(r,j-1)+AS1(r,j)*T2(r-1,j))./AP1(r,j);
? ? ? ? ? ?end
? ? ? ? ? ?for r=sizey:sizey%中部頂邊
? ? ? ? ? ? ? ?B(r,j)=s*dx*dy+ap0*Tmatrix(i-1,r,j)+TN*AN1(r,j);
? ? ? ? ? ? ? ?T2(r,j)=(B(r,j)+AE1(r,j)*T1(r,j+1)+AW1(r,j)*T2(r,j-1)+AS1(r,j)*T2(r-1,j))./AP1(r,j);
? ? ? ? ? ?end
? ? ? ?end
? ? ? ?for j=sizex:sizex%右側(cè)
? ? ? ? ? ?for r=1:1%右下角
? ? ? ? ? ? ? ?B(r,j)=s*dx*dy+ap0*Tmatrix(i-1,r,j)+TE*AE1(r,j)+TS*AS1(r,j);
? ? ? ? ? ? ? ?T2(r,j)=(B(r,j)+AW1(r,j)*T2(r,j-1)+AN1(r,j)*T1(r+1,j))./AP1(r,j);
? ? ? ? ? ?end
? ? ? ? ? ?for r=2:sizey-1%右側(cè)中部
? ? ? ? ? ? ? ?B(r,j)=s*dx*dy+ap0*Tmatrix(i-1,r,j)+TE*AE1(r,j);
? ? ? ? ? ? ? ?T2(r,j)=(B(r,j)+AW1(r,j)*T2(r,j-1)+AS1(r,j)*T2(r-1,j)+AN1(r,j)*T1(r+1,j))./AP1(r,j);
? ? ? ? ? ?end
? ? ? ? ? ?for r=sizey:sizey%右上角
? ? ? ? ? ? ? ?B(r,j)=s*dx*dy+ap0*Tmatrix(i-1,r,j)+TE*AE1(r,j)+TN*AN1(r,j);
? ? ? ? ? ? ? ?T2(r,j)=(B(r,j)+AS1(r,j)*T2(r-1,j)+AW1(r,j)*T2(r,j-1))./AP1(r,j);
? ? ? ? ? ?end
? ? ? ?end
? ? ? ?diff=abs(T2-T1);
? ? ? ?Err=max(abs(diff(:)));
? ? ? ?T3=T1;
? ? ? ?T1=T2;
? ?end
? ?for r=1:sizey
? ? ? ?for j=1:sizex
? ? ? ? ? ?Tmatrix(i,r,j)=T2(r,j);
? ? ? ?end
? ?end
? ?%把結(jié)果保存進(jìn)Tmatrix的這一時(shí)層，方便下一個(gè)時(shí)間步的迭代
? ?T4=Tmatrix(i,:,:);
? ?i
? ?%輸出一下i，看看算到第幾步了。
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
