传热学3-10
n=10 % x节点数
tm=480 % 时间节点
tfa=20; % 室内温度
np=10;
Bia=6*0.24/(n-1)/0.43; % 计算Bi及Fo准则数初始值
Bib=15*0.24/(n-1)/0.43;
Fo=0.43/1668/750*86400/tm/(0.24/(n-1))^2;
p=(2*Bib+2)*Fo % 显式差分格式的稳定性条件P<1
t=ones(1,n); % 输入墙体内部温度的初始值，各节点均设为1
tw=ones(1,tm+1);
tn=ones(1,tm+1);
if(p>1) % 迭代程序开始
disp('不稳定')
else
for s=1:tm+1
tfb=-3-6*cos((s-1)/tm*2*pi); % 室外温度
for k=2:n-1
t(1)=2*Fo*(t(2)+tfa*Bia)+(1-2*Fo-2*Bia*Fo)*t(1);
t(n)=2*Fo*(t(n-1)+tfb*Bib)+(1-2*Fo-2*Bib*Fo)*t(n-1);
t(k)=Fo*(t(k-1)+t(k+1))+(1-2*Fo)*t(k);
end
tn(s)=t(1);
tw(s)=t(n);
end % 迭代程序结束
Bia % 输出计算结果
Bib
Fo
for k=1:tm/np
tnn(1,k)=tn(1,np*k);
tww(1,k)=tw(1,np*k);
end
tnn
tww
w=0:tm; % 绘制温度变化图
plot(24*w/tm,tn(w+1),'-r',24*w/tm,tw(w+1)),grid
axis([0,24,-10,22])
set(gca,'xtick',[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24])
set(gca,'ytick',[-10,-8,-6,-4,-2,0,2,4,6,8,10,12,14,16,18,20,22,24]);
xlabel('时间'),ylabel('温度（℃）')
legend('内壁面温度','外壁面温度',0)
end
==============================================================================
clear;clc;
m=10
tm=480
t=ones(1,m);
tw=ones(1,tm+1);
tn=ones(1,tm+1);
tfa=20;
Bia=6*0.24/(m-1)/0.43;
Bib=15*0.24/(m-1)/0.43;
Fo=0.43/1668/750*86400/tm/(0.24/(m-1))^2;
p=(2*Bib+2)*Fo
for s=1:tm+1
tfb1=-9+(s-1)/tm*42;
tfb2=12-(s-1-tm/2)/tm*42;
tfbb(s)=tfb1
tfbc(s)=tfb2
for n=1:500
for k=2:m-1
t(1)=2*Fo*(t(2)+tfa*Bia)+(1-2*Fo-2*Bia*Fo)*t(1);
if s
else
t(m)=2*Fo*(t(m-1)+tfb2*Bib)+(1-2*Fo-2*Bib*Fo)*t(m-1);
end
t(k)=Fo*(t(k-1)+t(k+1))+(1-2*Fo)*t(k);
end
tn(s)=t(1);
tw(s)=t(m);
end
end
Bia
Bib
Fo
t
w=0:tm;
plot(24*w/tm,tn(w+1),'-r',24*w/tm,tw(w+1)),grid
axis([0,24,-2,22])
set(gca,'xtick',[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24])
set(gca,'ytick',[-2,0,2,4,6,8,10,12,14,16,18,20,22,24]);
xlabel('时间'),ylabel('温度（℃）')
legend('内壁面温度','外壁面温度',0)