二維非穩(wěn)態(tài)導(dǎo)熱的ADI迭代格式Matlab代碼——第一類邊界條件
%%感謝up主 @計(jì)算傳熱學(xué)大叔 的視頻講解?。?！
%% 線迭代基本思想：在每個(gè)時(shí)間步內(nèi)的每個(gè)迭代層中對矩形計(jì)算域從下往上逐行TDMA，
% 對于第n迭代層內(nèi)的行內(nèi)元素，Ten與Twn是未知量，與Tpn組成三對角矩陣的解
% 行下的Tsn已經(jīng)知道了，行上的元素Tnn不知道，所以用上個(gè)迭代層得到的Tnn-1去代替
% 參考1D隱式求解程序，每一行元素T的TDMA需要一個(gè)sizex維方陣存放系數(shù)
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維矩陣——對每個(gè)時(shí)間步，每一"行（x）"網(wǎng)格對應(yīng)一個(gè)對角陣，對角矩陣沿對角線連接。
%這樣做的好處是可以將整個(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ù)矩陣
%沒設(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
%% 開始迭代
for i=2:sizet
? ?T1=T0+zeros(sizey,sizex);%迭代初值
? ?T2=T1;
? ?ErrTol=1e-6;
? ?Err=1;
? ?while Err>ErrTol
? ? ? ?%隱式線迭代我想用東西線從南往北迭代，那么就必須從第一行往最后一行寫，使得Ts能一直迭代
? ? ? ?for j=1:1%第一行
? ? ? ? ? ?for r=1:1%左下角
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TW*AW1(j,r)+TS*AS1(j,r)+AN1(j,r).*T1(j+1,r);
? ? ? ? ? ?end
? ? ? ? ? ?for r=2:sizex-1%底部中部
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TS*AS1(j,r)+AN1(j,r)*T1(j+1,r);
? ? ? ? ? ?end
? ? ? ? ? ?for r=sizex:sizex%右下角
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TE*AE1(j,r)+TS*AS1(j,r)+AN1(j,r)*T1(j+1,r);
? ? ? ? ? ?end
? ? ? ? ? ?A1=zeros(sizex,sizex);
? ? ? ? ? ?b1=zeros(sizex,1);
? ? ? ? ? ?for r=1:sizex
? ? ? ? ? ? ? ?b1(r)=B(j,r);
? ? ? ? ? ? ? ?A1(r,r)=AP1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for r=1:sizex-1
? ? ? ? ? ? ? ?A1(r,r+1)=-AE1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for r=2:sizex
? ? ? ? ? ? ? ?A1(r,r-1)=-AW1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?%系數(shù)矩陣和常數(shù)向量確定完畢
? ? ? ? ? ?T2(j,:)=tdma(A1,b1);
? ? ? ?end
? ? ? ?for j=2:sizey-1
? ? ? ? ? ?for r=1:1%左邊
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TW*AW1(j,r)+AN1(j,r)*T1(j+1,r)+AS1(j,r)*T2(j-1,r);
? ? ? ? ? ? ? ?%千萬注意AS1(j,r)*T2(j-1,r),因?yàn)槭菑哪贤彼?，T2(j-1,r)已知
? ? ? ? ? ?end
? ? ? ? ? ?for r=2:sizex-1%主體
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+AN1(j,r)*T1(j+1,r)+AS1(j,r)*T2(j-1,r);
? ? ? ? ? ? ? ?%千萬注意AS1(j,r)*T2(j-1,r),因?yàn)槭菑哪贤彼悖琓2(j-1,r)已知
? ? ? ? ? ?end
? ? ? ? ? ?for r=sizex:sizex%右邊
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TE*AE1(j,r)+AS1(j,r)*T2(j-1,r)+AN1(j,r)*T1(j+1,r);
? ? ? ? ? ?end
? ? ? ? ? ?A2=zeros(sizex,sizex);
? ? ? ? ? ?b2=zeros(sizex,1);
? ? ? ? ? ?for r=1:sizex
? ? ? ? ? ? ? ?b2(r)=B(j,r);
? ? ? ? ? ? ? ?A2(r,r)=AP1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for r=1:sizex-1
? ? ? ? ? ? ? ?A2(r,r+1)=-AE1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for r=2:sizex
? ? ? ? ? ? ? ?A2(r,r-1)=-AW1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?T2(j,:)=tdma(A2,b2);
? ? ? ?end
? ? ? ?for j=sizey:sizey
? ? ? ? ? ?for r=1:1%左上角
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TW*AW1(j,r)+TN*AN1(j,r)+AS1(j,r)*T2(j-1,r);
? ? ? ? ? ?end
? ? ? ? ? ?for r=2:sizex-1%頂部中
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TN*AN1(j,r)+AS1(j,r)*T2(j-1,r);
? ? ? ? ? ?end
? ? ? ? ? ?for r=sizex:sizex%右上角
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TE*AE1(j,r)+TN*AN1(j,r)+AS1(j,r)*T2(j-1,r);
? ? ? ? ? ?end
? ? ? ? ? ?A3=zeros(sizex,sizex);
? ? ? ? ? ?b3=zeros(sizex,1);
? ? ? ? ? ?for r=1:sizex
? ? ? ? ? ? ? ?b3(r)=B(j,r);
? ? ? ? ? ? ? ?A3(r,r)=AP1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for r=1:sizex-1
? ? ? ? ? ? ? ?A3(r,r+1)=-AE1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for r=2:sizex
? ? ? ? ? ? ? ?A3(r,r-1)=-AW1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?T2(j,:)=tdma(A3,b3);
? ? ? ?end
? ? ? ?%% 到此新溫度場T2就得到了，ADI迭代就是在這里再補(bǔ)一次y方向的線迭代得到溫度場T3，最終比較T3和T1的大小
? ? ? ?for r=1:1%從左下往右上迭代
? ? ? ? ? ?for j=1:1%左下角
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TW*AW1(j,r)+TS*AS1(j,r)+AE1(j,r)*T2(j,r+1);
? ? ? ? ? ?end
? ? ? ? ? ?for j=2:sizey-1%左中
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TW*AW1(j,r)+AE1(j,r)*T2(j,r+1);
? ? ? ? ? ?end
? ? ? ? ? ?for j=sizey:sizey%左上
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TW*AW1(j,r)+TN*AN1(j,r)+AE1(j,r)*T2(j,r+1);
? ? ? ? ? ?end
? ? ? ? ? ?T3=T2;
? ? ? ? ? ?A4=zeros(sizey,sizey);
? ? ? ? ? ?b4=zeros(sizey,1);
? ? ? ? ? ?for j=1:sizey
? ? ? ? ? ? ? ?b4(j)=B(j,r);
? ? ? ? ? ? ? ?A4(j,j)=AP1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for j=1:sizey-1
? ? ? ? ? ? ? ?A4(j,j+1)=-AN1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for j=2:sizey
? ? ? ? ? ? ? ?A4(j,j-1)=-AS1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?T3(:,r)=tdma(A4,b4);
? ? ? ?end
? ? ? ?for r=2:sizex-1
? ? ? ? ? ?for j=1:1%中下
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TS*AS1(j,r)+AE1(j,r)*T2(j,r+1)+AW1(j,r)*T2(j,r-1);
? ? ? ? ? ?end
? ? ? ? ? ?for j=2:sizey-1%主體
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+AE1(j,r)*T2(j,r+1)+AW1(j,r)*T2(j,r-1);
? ? ? ? ? ?end
? ? ? ? ? ?for j=sizey:sizey%中上
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TN*AN1(j,r)+AE1(j,r)*T2(j,r+1)+AW1(j,r)*T2(j,r-1);
? ? ? ? ? ?end
? ? ? ? ? ?A5=zeros(sizey,sizey);
? ? ? ? ? ?b5=zeros(sizey,1);
? ? ? ? ? ?for j=1:sizey
? ? ? ? ? ? ? ?b5(j)=B(j,r);
? ? ? ? ? ? ? ?A5(j,j)=AP1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for j=1:sizey-1
? ? ? ? ? ? ? ?A5(j,j+1)=-AN1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for j=2:sizey
? ? ? ? ? ? ? ?A5(j,j-1)=-AS1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?T3(:,r)=tdma(A5,b5);
? ? ? ?end
? ? ? ?for r=sizex:sizex
? ? ? ? ? ?for j=1:1%右下
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TE*AE1(j,r)+TS*AS1(j,r)+AW1(j,r)*T2(j,r-1);
? ? ? ? ? ?end
? ? ? ? ? ?for j=2:sizey-1%右邊
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TE*AE1(j,r)+AW1(j,r)*T2(j,r-1);
? ? ? ? ? ?end
? ? ? ? ? ?for j=sizey:sizey%右上
? ? ? ? ? ? ? ?B(j,r)=s*dx*dy+ap0*Tmatrix(i-1,j,r)+TE*AE1(j,r)+TN*AN1(j,r)+AW1(j,r)*T2(j,r-1);
? ? ? ? ? ?end
? ? ? ? ? ?A6=zeros(sizey,sizey);
? ? ? ? ? ?b6=zeros(sizey,1);
? ? ? ? ? ?for j=1:sizey
? ? ? ? ? ? ? ?b6(j)=B(j,r);
? ? ? ? ? ? ? ?A6(j,j)=AP1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for j=1:sizey-1
? ? ? ? ? ? ? ?A6(j,j+1)=-AN1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?for j=2:sizey
? ? ? ? ? ? ? ?A6(j,j-1)=-AS1(j,r);
? ? ? ? ? ?end
? ? ? ? ? ?T3(:,r)=tdma(A6,b6);
? ? ? ?end
? ? ? ?diff=abs(T3-T1);
? ? ? ?Err=max(abs(diff(:)));
? ? ? ?T3=T1;
? ? ? ?T1=T2;
? ?end
? ?for r=1:sizey
? ? ? ?for j=1:sizex
? ? ? ? ? ?Tmatrix(i,r,j)=T3(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
function x=tdma(A,b)
? ?n=length(b);
? ?alpha=zeros(n,1);
? ?beta=zeros(n,1);
? ?% 前向消去
? ?alpha(2)=-A(1,2)/A(1,1);
? ?beta(2)=b(1)/A(1,1);
? ?for i=2:n-1
? ? ? ?alpha(i+1)=-A(i,i+1)/(A(i,i)+A(i,i-1)*alpha(i));
? ? ? ?beta(i+1)=(b(i)-A(i,i-1)*beta(i))/(A(i,i)+A(i,i-1)*alpha(i));
? ?end
? ?% 回代求解
? ?x=zeros(n,1);
? ?x(n)=(b(n)-A(n,n-1)*beta(n))/(A(n,n)+A(n,n-1)*alpha(n));
? ?for i=n-1:-1:1
? ? ? ?x(i)=alpha(i+1)*x(i+1)+beta(i+1);
? ?end
end
%這是標(biāo)準(zhǔn)的TDMA解方程程序，調(diào)用格式為T=tdma(A,b)
