二维时间隐式格式高斯-赛德尔(G-S)点迭代MATLAB程序

最新推荐文章于 2022-04-19 20:24:23 发布

原创最新推荐文章于 2022-04-19 20:24:23 发布 · 2.5k 阅读

14 ·

CC 4.0 BY-SA版权

文章标签：

#matlab

数值传热学专栏收录该内容

7 篇文章

订阅专栏

本文提供了一个使用MATLAB实现的二维时间隐式格式高斯-赛德尔(G-S)点迭代法程序示例。该程序用于解决数值传热学问题，通过设置边界条件和网格系数，进行迭代求解直至收敛。适用于传热学研究及教学。

二维时间隐式格式高斯-赛德尔(G-S)点迭代MATLAB程序

代码改编自计算传热学–15-高斯-赛德尔(G-S)点迭代-数值传热学
下面是代码

%GSiteration implicit time 2d
clc
clear
%%
LengthX = 3;
LengthY = 2;
Tleft = 3.5;
Tright = 5;
Tbottom = 2;
Tup = 4;
den = 100;
s = 10;
k=10;
cv = 1000;
%%
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);

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;
    for iter = 1:itermax
        for i=M+1:-1:2
            for j=2:N+1
                T1(i,j) = (ae1(i,j)*Tn0(i,j+1)+aw1(i,j)*T1(i,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))/ap1(i,j);
            end
        end
        resimax = max(max(abs(T1-Tn0)));
        if resimax<maxresi
            break;
        end
        Tn0 = T1;
    end
    resiall(h) = resimax;
    iterall(h) = iter;
    T0=T1;
    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)+2,2);
end
plot(dt:dt:dt*maxstep,T11(2:maxstep+1));
hold on