非常感谢“计算传热学大叔”，大家了解更多，请移步前期文章：https://blog.csdn.net/weixin_37928884/article/details/127709215
第一类边界条件
clc
clear
close all %直接在此修改参数
length = 0.135; %长度
Tb = 930; %初始温度
TL = 1150; %边界温度1
TR = 1150; %边界温度2
den = 7840; %控制体密度
C = 465; %控制体比热
k = 28.5; %控制体导热系数
s = 0; %源项
h = 375; %换热系数
dt = 1; %步长
n = 10; %网格数
steps =3000; %步数
dx = length / n;
ae0 = zeros(1,n);
aw0 = zeros(1,n);
ap0 = zeros(1,n);
ap1 = zeros(1,n);
T0 = Tb * ones(steps+1 ,n+2);
T0(:,1) = TL * ones(1,steps + 1).';
T0(:,n + 2) = TR * ones(1,steps + 1).';
b = s * dx;
for x = 1:steps
for i = 2:n
ae0(1,i) = k / dx;
aw0(1,i) = k / dx;
ap0(1,i) = den * C * dx / dt - ae0(1,i) - aw0(1,i);
ap1(1,i) = ae0(1,i) + aw0(1,i) + ap0(1,i);
end
for i = 1
ae0(1,i) = k /( dx / 2);
aw0(1,i) = k / dx;
ap0(1,i) = den * C * dx / dt - ae0(1,i) - aw0(1,i);
ap1(1,i) = ae0(1,i) + aw0(1,i) + ap0(1,i);
end
for i = n
ae0(1,i) = k / dx;
aw0(1,i) = k /( dx / 2);
ap0(1,i) = den * C * dx / dt - ae0(1,i) - aw0(1,i);
ap1(1,i) = ae0(1,i) + aw0(1,i) + ap0(1,i);
end
end
for x = 2:steps+1
for i = 2 : n + 1
T0(x , i) = ( ae0(1,i-1) * T0(x-1,i+1) + aw0(1,i-1) * T0 (x-1,i-1) + ap0(1,i-1)* T0(x-1,i) + b)/ap1(1,i-1);
end
end
xlist = length/(2 * n) : length/n : length - length/(2 * n); %坐标（不包括边界）
Tend = T0(steps + 1,2 : n +1); %结束时温度
ylist = T0(steps + 1 ,:);
result =[0 xlist length ylist]; %坐标和结束温度
subplot(1,2,1);
plot(xlist,Tend)
xlabel('坐标/m');
ylabel('温度/℃');
title(['末态温度分布/网格数:',num2str(n)]);
grid on;
subplot(1,2,2);
for i = 1 : n + 1
x = 1 : steps + 1;
time = x * dt/3600;
y = T0(:,i);
plot(time,y);
title(['各点温度与时间关系/步数：',num2str(steps),'/步长:',num2str(dt),'s']);
xlabel('时间/h');
ylabel('温度/℃');
hold on
end
第三类边界条件
clc
clear
close all %直接在此修改参数
length = 0.135; %长度
Tb = 930; %初始温度
TL = 1150; %边界温度1
TR = 1150; %边界温度2
den = 7840; %控制体密度
C = 465; %控制体比热
k = 28.5; %控制体导热系数
s = 0; %源项
h = 375; %换热系数
dt = 1; %步长
n = 10; %网格数
steps =3000; %步数
dx = length / n;
ae0 = zeros(1,n);
aw0 = zeros(1,n);
ap0 = zeros(1,n);
ap1 = zeros(1,n);
T0 = Tb * ones(steps+1 ,n+2);
T0(:,1) = TL * ones(1,steps + 1).';
T0(:,n + 2) = TR * ones(1,steps + 1).';
b = s * dx;
for x = 1:steps
for i = 2:n
ae0(1,i) = k / dx;
aw0(1,i) = k / dx;
ap0(1,i) = den * C * dx / dt - ae0(1,i) - aw0(1,i);
ap1(1,i) = ae0(1,i) + aw0(1,i) + ap0(1,i);
end
for i = 1
ae0(1,i) = 1/(((dx / 2)/k) + 1/h) ;%k /( dx / 2);
aw0(1,i) = k / dx;
ap0(1,i) = den * C * dx / dt - ae0(1,i) - aw0(1,i);
ap1(1,i) = ae0(1,i) + aw0(1,i) + ap0(1,i);
end
for i = n
ae0(1,i) = k / dx;
aw0(1,i) = 1/(((dx / 2)/k) + 1/h);%k /( dx / 2);
ap0(1,i) = den * C * dx / dt - ae0(1,i) - aw0(1,i);
ap1(1,i) = ae0(1,i) + aw0(1,i) + ap0(1,i);
end
end
for x = 2:steps+1
for i = 2 : n + 1
T0(x , i) = ( ae0(1,i-1) * T0(x-1,i+1) + aw0(1,i-1) * T0 (x-1,i-1) + ap0(1,i-1)* T0(x-1,i) + b)/ap1(1,i-1);
end
end
xlist = length/(2 * n) : length/n : length - length/(2 * n); %坐标（不包括边界）
Tend = T0(steps + 1,2 : n +1); %结束时温度
ylist = T0(steps + 1 ,:);
result =[0 xlist length ylist]; %坐标和结束温度
subplot(1,2,1);
plot(xlist,Tend)
xlabel('坐标/m');
ylabel('温度/℃');
title(['末态温度分布/网格数:',num2str(n)]);
grid on;
subplot(1,2,2);
for i = 1 : n + 1
x = 1 : steps + 1;
time = x * dt/3600;
y = T0(:,i);
plot(time,y);
title(['各点温度与时间关系/步数：',num2str(steps),'/步长:',num2str(dt),'s']);
xlabel('时间/h');
ylabel('温度/℃');
hold on
end
