MATLAB一维瞬态导热时间半隐格式

代码改变自计算传热学–8-时间半隐格式离散方程和统一的时间显式/隐式/半隐格式C程序-数值传热学_哔哩哔哩 (゜-゜)つロ 干杯~-bilibili

% implicitExplicit time
%%
clc
clear
Llength = 3;
Tleft = 3;
Tright = 5;
density = 100;
source = 10;
c = 1000;
deltaT = 200;
k=10;
N = 5;
maxstep = 2000;
%%
ae0 = zeros(N+2,1);
ae1 = zeros(N+2,1);
aw0 = zeros(N+2,1);
aw1 = zeros(N+2,1);
ap0 = zeros(N+2,1);
ap1 = zeros(N+2,1);
b = zeros(N+2,1);
T1 = zeros(N+2,1);
a = zeros(N+2,1);
b1 = zeros(N+2,1);
c1 = zeros(N+2,1);
d = zeros(N+2,1);
Tall = zeros(N+2,maxstep+1);
%%
%内部网格计算
dx = Llength/N;
ae0(3:N,1) = k/(2*dx);
aw0(3:N,1) = k/(2*dx);
ae1(3:N,1) = k/(2*dx);
aw1(3:N,1) = k/(2*dx);
ap0(3:N,1) = density*c*dx/deltaT - k/(2*dx) - k/(2*dx);
ap1(3:N,1) = ae0(3:N,1)+aw0(3:N,1)+aw1(3:N,1)+ae1(3:N,1)+ap0(3:N,1);
b(3:N,1) = source*dx;
%边界网格计算
%第一个
i=2;
ae0(i,1) = k/(2*dx);
aw0(i,1) = k/(2*dx/2);
ae1(i,1) = k/(2*dx);
aw1(i,1) = k/(2*dx/2);
ap0(i,1) = density*c*dx/deltaT - k/(2*dx/2) - k/(2*dx);
ap1(i,1) = ae0(i,1)+aw0(i,1)+aw1(i,1)+ae1(i,1)+ap0(i,1);
b(i,1) = source*dx;
%最后一个
i =N+1;
ae0(i,1) = k/(2*dx/2);
aw0(i,1) = k/(2*dx);
ae1(i,1) = k/(2*dx/2);
aw1(i,1) = k/(2*dx);
ap0(i,1) = density*c*dx/deltaT - k/(2*dx) - k/(2*dx/2);
ap1(i,1) = ae0(i,1)+aw0(i,1)+aw1(i,1)+ae1(i,1)+ap0(i,1);
b(i,1) = source*dx;
%初始条件
T0(2:N+1,1) = Tleft;
T0(1,1)=Tleft;
T0(N+2,1)=Tright;
T1(1,1) = Tleft;
T1(N+2,1)=Tright;
Tall(:,1) = T0;
%abcd初始化
for i =1:maxstep
    g=2;
    a(g,1) = ap1(g,1);
    b1(g,1) = -aw1(g,1);
    c1(g,1) = -ae1(g,1);
    d(g,1) = ap0(g,1)*T0(g,1) + aw1(g,1)*T1(g-1,1)+b(g,1)+ae0(g,1)*T0(g+1,1)+aw0(g,1)*T0(g-1,1);
    
    a(3:N,1) = ap1(3:N,1);
    b1(3:N,1) = -aw1(3:N,1);
    c1(3:N,1) = -ae1(3:N,1);
    d(3:N,1) = ap0(3:N,1).*T0(3:N,1) +b(3:N,1)+ae0(3:N,1).*T0(3+1:N+1,1)+aw0(3:N,1).*T0(3-1:N-1,1);
    %第一个
    
      g=N+1;
    a(g,1) = ap1(g,1);
    b1(g,1) = -aw1(g,1);
    c1(g,1) = -ae1(g,1);
    d(g,1) = ae1(g,1)*T1(g+1,1)+ap0(g,1)*T0(g,1) +b(g,1)+ae0(g,1)*T0(g+1,1)+aw0(g,1)*T0(g-1,1);
    %计算T1时刻
    [T1]=TMDA(a,b1,c1,d,N,Tleft,Tright);
    
    T0 = T1;
    Tall(:,i+1) = T0;
end
Tall = Tall';
plot(Tall(:,2:N+1));
lgd =legend('T1','T2','T3','T4','T5');
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