二维瞬态导热时间隐式格式ADI迭代MATLAB程序

二维瞬态导热时间隐式格式ADI迭代MATLAB程序

代码改编自计算传热学–17-ADI迭代-数值传热学
计算传热学–17-ADI迭代的源程序
下面是代码

%ADI implicit time 2d
clc
clear
%%
LengthX = 3;
LengthY = 2;
Tleft = 3.5;
Tright = 5;
Tbottom = 2;
Tup = 4;
den = 100;
cv = 1000;
s = 0;
k=10;
%%
dt = 200;
N = 15;%X
M = 15;%Y
maxstep = 200;
itermax = 1000;
maxresi = 1e-7;
%%
dx =LengthX/N;
dy =LengthY/M;
ae1 = zeros(M+2,N+2);
aw1 = zeros(M+2,N+2);
an1 = zeros(M+2,N+2);
as1 = zeros(M+2,N+2);
ap0 = zeros(M+2,N+2);
ap1 = zeros(M+2,N+2);
b= zeros(M+2,N+2);
T0 =zeros(M+2,N+2);
T1 = zeros(M+2,N+2);

a = zeros(N+2,1);
b1 = zeros(N+2,1);
c1 = zeros(N+2,1);
d = zeros(N+2,1);

resi = zeros(M+2,N+2);
Tall = zeros(M+2,maxstep+1);
resiall = zeros(maxstep,1);
iterall = zeros(maxstep,1);
%% 边界条件
T0(2:M+1,2:N+1) = 3;
T0(1:M+2,1) = Tleft;
T0(1:M+2,N+2) = Tright;
T0(1,1:N+2) = Tup;
T0(M+2,1:N+2) = Tbottom;

T1(1:M+2,1) = Tleft;
T1(1:M+2,N+2) = Tright;
T1(1,1:N+2) = Tup;
T1(M+2,1:N+2) = Tbottom;
%% 网格系数计算
%内部网格计算
ae1(3:M,3:N) = k*dy/dx;
aw1(3:M,3:N) = k*dy/dx;
an1(3:M,3:N) = k*dx/dy;
as1(3:M,3:N) = k*dx/dy;
ap0(3:M,3:N) = den*cv*dx*dy/dt;
ap1(3:M,3:N) = ap0(3:M,3:N)+ae1(3:M,3:N)+aw1(3:M,3:N)+an1(3:M,3:N)+as1(3:M,3:N);
b(3:M,3:N) = s*dx*dy;
%边界网格计算
%第一个
i=2;
for j=3:N
    ae1(i,j) = k*dy/dx;
    aw1(i,j) = k*dy/dx;
    an1(i,j) = k*dx/(dy/2);
    as1(i,j) = k*dx/dy;
    ap0(i,j) = den*cv*dx*dy/dt;
    ap1(i,j) = ap0(i,j)+ae1(i,j)+aw1(i,j)+an1(i,j)+as1(i,j);
    b(i,j) = s*dx*dy;
end

i=M+1;
for j=3:N
    ae1(i,j) = k*dy/dx;
    aw1(i,j) = k*dy/dx;
    an1(i,j) = k*dx/dy;
    as1(i,j) = k*dx/(dy/2);
    ap0(i,j) = den*cv*dx*dy/dt;
    ap1(i,j) = ap0(i,j)+ae1(i,j)+aw1(i,j)+an1(i,j)+as1(i,j);
    b(i,j) = s*dx*dy;
end
j=2;
for i=3:M
    ae1(i,j) = k*dy/dx;
    aw1(i,j) = k*dy/(dx/2);
    an1(i,j) = k*dx/dy;
    as1(i,j) = k*dx/dy;
    ap0(i,j) = den*cv*dx*dy/dt;
    ap1(i,j) = ap0(i,j)+ae1(i,j)+aw1(i,j)+an1(i,j)+as1(i,j);
    b(i,j) = s*dx*dy;
end
j=N+1;
for i=3:M
    ae1(i,j) = k*dy/(dx/2);
    aw1(i,j) = k*dy/dx;
    an1(i,j) = k*dx/dy;
    as1(i,j) = k*dx/dy;
    ap0(i,j) = den*cv*dx*dy/dt;
    ap1(i,j) = ap0(i,j)+ae1(i,j)+aw1(i,j)+an1(i,j)+as1(i,j);
    b(i,j) = s*dx*dy;
end

j=2;
i=2;
ae1(i,j) = k*dy/dx;
aw1(i,j) = k*dy/(dx/2);
an1(i,j) = k*dx/(dy/2);
as1(i,j) = k*dx/dy;
ap0(i,j) = den*cv*dx*dy/dt;
ap1(i,j) = ap0(i,j)+ae1(i,j)+aw1(i,j)+an1(i,j)+as1(i,j);
b(i,j) = s*dx*dy;

j=2;
i=M+1;
ae1(i,j) = k*dy/dx;
aw1(i,j) = k*dy/(dx/2);
an1(i,j) = k*dx/dy;
as1(i,j) = k*dx/(dy/2);
ap0(i,j) = den*cv*dx*dy/dt;
ap1(i,j) = ap0(i,j)+ae1(i,j)+aw1(i,j)+an1(i,j)+as1(i,j);
b(i,j) = s*dx*dy;

i=2;
j=N+1;
ae1(i,j) = k*dy/(dx/2);
aw1(i,j) = k*dy/dx;
an1(i,j) = k*dx/(dy/2);
as1(i,j) = k*dx/dy;
ap0(i,j) = den*cv*dx*dy/dt;
ap1(i,j) = ap0(i,j)+ae1(i,j)+aw1(i,j)+an1(i,j)+as1(i,j);
b(i,j) = s*dx*dy;

i=M+1;
j=N+1;
ae1(i,j) = k*dy/(dx/2);
aw1(i,j) = k*dy/dx;
an1(i,j) = k*dx/dy;
as1(i,j) = k*dx/(dy/2);
ap0(i,j) = den*cv*dx*dy/dt;
ap1(i,j) = ap0(i,j)+ae1(i,j)+aw1(i,j)+an1(i,j)+as1(i,j);
b(i,j) = s*dx*dy;
%%
Tall(1:M+2,1:N+2) = T0;
%% iteration

for h = 1:maxstep
    Tn0 = T0;
    %x方向
    for iter=1:itermax
        for i =M+1:-1:2
            %             i=i;
            j=2;
            a(j,1) = ap1(i,j);
            b1(j,1) = -aw1(i,j);
            c1(j,1) = -ae1(i,j);
            d(j,1) = an1(i,j)*Tn0(i-1,j)+as1(i,j)*T1(i+1,j)+ap0(i,j)*T0(i,j) + aw1(i,j)*T1(i,j-1)+b(i,j);
            for j=3:N
                a(j,1) = ap1(i,j);
                b1(j,1) = -aw1(i,j);
                c1(j,1) = -ae1(i,j);
                d(j,1) = an1(i,j)*Tn0(i-1,j)+as1(i,j)*T1(i+1,j)+ap0(i,j)*T0(i,j) + b(i,j);
            end
            j=N+1;
            a(j,1) = ap1(i,j);
            b1(j,1) = -aw1(i,j);
            c1(j,1) = -ae1(i,j);
            d(j,1) = an1(i,j)*Tn0(i-1,j)+as1(i,j)*T1(i+1,j)+ap0(i,j)*T0(i,j) + ae1(i,j)*T1(i,j+1)+b(i,j);
            [x]=TMDA(a,b1,c1,d,N,Tleft,Tright);
            T1(i,:)=x;
            
        end
        %Y方向
        Tn0 = T1;
        for j =2:N+1
            %             i=i;
            i=N+1;
            a(2,1) = ap1(i,j);
            %         b1(i,1) = -as1(i,j);
            c1(2,1) = -an1(i,j);
            d(2,1) = ae1(i,j)*Tn0(i,j+1)+aw1(i,j)*T1(i,j-1)+ap0(i,j)*T0(i,j) + as1(i,j)*T1(i+1,j)+b(i,j);
            for i=M:-1:3
                t = M-i+3;
                a(t,1) = ap1(i,j);
                b1(t,1) = -as1(i,j);
                c1(t,1) = -an1(i,j);
                d(t,1) = ae1(i,j)*Tn0(i,j+1)+aw1(i,j)*T1(i,j-1)+ap0(i,j)*T0(i,j) +b(i,j);
            end
            i=2;
            a(N+1,1) = ap1(i,j);
            b1(N+1,1) = -as1(i,j);
            %         c1(i,1) = -an1(i,j);
            d(N+1,1) = ae1(i,j)*Tn0(i,j+1)+aw1(i,j)*T1(i,j-1)+ap0(i,j)*T0(i,j) + an1(i,j)*T1(i-1,j)+b(i,j);
            [x]=TMDA(a,b1,c1,d,M,Tbottom,Tup);
            T1(:,j)=flipud(x);
            
        end
        resimax = max(max(abs(T1-Tn0)));
        if resimax<maxresi
            break;
        end
        Tn0 = T1;
    end
    T0=T1;
    resiall(h) = resimax;
    iterall(h) = iter;
    
    Tall((h)*(M+2)+1:(h+1)*(M+2),1:N+2) = T0;
end




% plot(linspace(0,LengthX,N),Tall(maxstep*(M+2)+2,2:N+1));
T11=zeros(maxstep,1);
for i =1:maxstep+1
    T11(i) = Tall((i-1)*(M+2)+6,2);
end
plot(dt:dt:dt*maxstep,T11(2:maxstep+1));
hold on


function [x]=TMDA(a,b1,c1,d,N,Tleft,Tright)
p = zeros(N+2,1);
q = zeros(N+2,1);
x = zeros(N+2,1);
i=2;
p(i,1) = -c1(i)/a(i);
q(i) = d(i)/a(i);
%中间
for i=3:N
    p(i,1) = -c1(i,1)/(a(i,1)+b1(i,1)*p(i-1,1));
    q(i,1) = (d(i,1)-b1(i,1)*q(i-1,1))/(a(i,1)+b1(i,1)*p(i-1,1));
end
%最后一个
i=N+1;
x(i,1) = (d(i,1)-b1(i,1)*q(i-1,1))/(a(i,1)+b1(i,1)*p(i-1,1));
for i=2:N
    x(N+2-i,1) = p(N+2-i,1)*x(N+3-i,1)+q(N+2-i,1);
end
x(1,1) = Tleft;
x(N+2,1)=Tright;
end