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

代码改编自计算传热学–6-一维瞬态常系数导热离散方程时间隐式格式

%%
clc
clear
Llength = 3;
Tleft = 3;
Tright = 5;
density = 100;
source = 10;
c = 1000;
deltaT = 200;
k=10;
N = 5;
maxstep = 1000;
%%
ae1 = 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;
ae1(3:N,1) = k/dx;
aw1(3:N,1) = k/dx;
ap0(3:N,1) = density*c*dx/deltaT;
ap1(3:N,1) = density*c*dx/deltaT+aw1(3:N,1)+ae1(3:N,1);
b(3:N,1) = source*dx;
%边界网格计算
%第一个
i=2;
ae1(i,1) = k/dx;
aw1(i,1) = k/(dx/2);
ap0(i,1) = density*c*dx/deltaT;
ap1(i,1) = density*c*dx/deltaT+aw1(i,1)+ae1(i,1);
b(i,1) = source*dx;
%最后一个
i=N+1;
ae1(i,1) = k/(dx/2);
aw1(i,1) = k/dx;
ap0(i,1) = density*c*dx/deltaT;
ap1(i,1) = ap0(i,1)+aw1(i,1)+ae1(i,1);
b(i,1) = source*dx;
%初始条件
T0(2:N+1,1) = Llength;
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);
    
    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);
      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);
    
    %计算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');
figure
plot(linspace(0,3,N+2),Tall(end,1:end));

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