二维非稳态导热的ADI迭代格式Matlab代码——第一类边界条件
%%感谢up主 @计算传热学大叔 的视频讲解！！！
%% 线迭代基本思想：在每个时间步内的每个迭代层中对矩形计算域从下往上逐行TDMA，
% 对于第n迭代层内的行内元素，Ten与Twn是未知量，与Tpn组成三对角矩阵的解
% 行下的Tsn已经知道了，行上的元素Tnn不知道，所以用上个迭代层得到的Tnn-1去代替
% 参考1D隐式求解程序，每一行元素T的TDMA需要一个sizex维方阵存放系数
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那样投机取巧——广义源项b必须是一个列向量才能求解方程，原有的三维数组/暴力求解方法失效
%可以试着降维成2维矩阵——对每个时间步，每一"行（x）"网格对应一个对角阵，对角矩阵沿对角线连接。
%这样做的好处是可以将整个计算域内的所有网格温度一起算，特别适合迭代求解。
Tmatrix=zeros(sizet,sizey,sizex);%时间，行，列
for i=1:1
Tmatrix(i,:,:)=T0;
end
%% 计算系数矩阵和源项矩阵
AP1=zeros(sizey,sizex);%ap1系数矩阵
AE1=zeros(sizey,sizex);%ae1系数矩阵
AW1=zeros(sizey,sizex);%aw1系数矩阵
AN1=zeros(sizey,sizex);%an1系数矩阵
AS1=zeros(sizey,sizex);%as1系数矩阵
%没设ap0，直接在b向量里计算就行
B=zeros(sizey,sizex);%广义源项矩阵
for j=1:1%按列去分
for r=1:1 %矩形区域左下顶点，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;
%广义源项需要用到上一时刻的温度，所以放到时间步内去算
end
for r=2:sizey-1 %矩形区域左边，dxw=0.5dx
%不用算ap0了，因为都一样
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 %矩形区域左上顶点，dxw=0.5dx,dyn=0.5dy
%不用算ap0了，因为都一样
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 %矩形区域中间下边，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 %矩形中间区域中间，主体部分
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 %矩形区域中间上边，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 %矩形区域右下顶点，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 %矩形区域右边，dxe=0.5dx
%不用算ap0了，因为都一样
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 %矩形区域右上顶点，dxw=0.5dx,dyn=0.5dy
%不用算ap0了，因为都一样
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
%系数矩阵和常数向量确定完毕
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),因为是从南往北算，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),因为是从南往北算，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)+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迭代就是在这里再补一次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
%把结果保存进Tmatrix的这一时层，方便下一个时间步的迭代
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,:)
% 创建网格坐标
[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
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
%这是标准的TDMA解方程程序，调用格式为T=tdma(A,b)
