节点 i 温度 t(℃)
1 250.00
2 232.01
3 213.59
4 194.71
5 175.35
6 155.46
7 134.98
8 113.87
9 92.06
10 69.47
11 46.00
热流密度 q= -4.47E+04 w/m2
（6） 温度分布图
10.一砖墙厚 240mm，内、外表面的表面传热系数分别为 6.0 w/(m2·k)和 15 w/(m2·k)，墙 3 体材料的导热系数λ=0.43 w/(m·k)，密度ρ=1668kg/m ，比热容 c=0.75kJ/(kg·K)，室内空 气温度保持不变为 20℃，室外空气温度周期性变化，中午 12 点温度最高为 3℃，晚上 12 点温度最低为-9℃，试用数值计算方法确定内、外墙壁面温度在一天内的变化。 解： （1） 建立模型 本题属于常物性，无热源的一维非稳态导热过程 将平壁沿厚度方向（x 方向）划分为 N 个均匀相等的间距，时间从τ=0 开始， 按Δτ分割为 k 段，令 i 表示内节点位置，k 表示 kΔτ时刻，节点布置如图所 示
开始
输入 n, np, tm, Bia, Bib, Fo
n=10, tm=480, tfa=20, Bia=6*0.24/(n-1)/0.43, Bib=15*0.24/(n-1)/0.43; Fo=0.43/1668/750*86400/tm/(0.24/(n-1))^2
No 打印“不稳定”
Yes k=k+1 k=1
1
2
3
4
N
N+1
本题给出 N=10。 （2） 通过热流法建立离散方程 a) 内节点离散方程
L
P
R
对节点 P(i)所代表的微元体，在 x 方向上与节点 P 相邻的节点分别为 L(i-1)和 R(i+1)。 由于节点之间的间距很小， 可以认为相邻节点间的温度分布是线性的。 节点 P 所代表的网格单元与它周围各网格单元之间的导热量可根据傅里叶定 律直接写为:
k ha ( t k f − t1 ) − λ k t1k − t2 t k +1 − t1k ∆x = ρc 1 ∆x ∆τ 2
整理上式，可得
ha ∆x k k 1 ρ c∆x 2 k +1 k t −t + ( t f − t1 ) = 2 λ∆τ ( t1 − t1 ) λ
k 2 k 1
上式中
τ k x i i=n+1
0 i=1
本题给出 N=10。 （2） 通过热流法建立离散方程（采用显式差分格式） a) 内节点离散方程 导热微分方程式为
∂t ∂ 2t =a 2 ∂τ ∂x
温度对时间的一阶导数，采用向前差分，则
(1)
tik +1 − tik ∂t = ∆τ ∂τ i ,k
外壁面温度(℃) -1.39 -1.86 -2.14 -2.30 -2.36 -2.32 -2.20 -2.01 -1.76 -1.47 -1.13 -0.76 -0.37 0.04 0.46 0.89 1.31 1.72 2.11 2.49 2.83 3.15 3.44 3.68
时间 12:30 13:00 13:30 14:00 14:30 15:00 15:30 16:00 16:30 17:00 17:30 18:00 18:30 19:00 19:30 20:00 20:30 21:00 21:30 22:00 22:30 23:00 23:30 0:00
Φ LP = λLP
其中
ti −1 − ti ×1 ∆x
Φ RP = λRP
ti +1 − ti ×1 ∆x
λRP = 43 + 0.08 × λLP
(ti + ti +1 ) 2 (t + t ) = 43 + 0.08 × i −1 i 2
对节点 P(i)所代表的微元写热平衡式，即可得节点 P(i)温度的离散方程
tfb=-3-6*cos((s-1)/tm*2*pi)
控制打印 NP No No k>tm
Yes
打印 tn, tw Yes
打印 Bia, Bib, Fo
停机
（4） 编写程序 本题使用 matlab 软件，所编写的程序如下： clear;clc; 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
温度对 x 的二阶导数，采用中心差分表达式应为
(2)
∂ 2t tik−1 − 2tik + tik+1 = 2 ∆x 2 ∂x i ,k tik +1 − tik t k − 2tik + tik+1 = a i −1 ∆τ ∆x 2
将上式移项整理得
(3)
将式(3)和式(2)代入式(1)，就得到内节点(i, k)的节点离散方程
Φ LP + Φ RP = 0
λLP ti =
b)
ti −1 − ti t −t ×1 + λRP i +1 i × 1 = 0 ∆x ∆x λLP ti −1 + λLP ti +1 λLP + λRP
边界节点离散方程 由于本题的壁面温度属于第一类边界条件，因此
t1 = 250 tn +1 = 46
（5） 计算结果 准则数初始值：Bia = 0.3721，Bib = 0.9302，Fo = 0.0870
时间 内壁面温度(℃) 0:30 8.22 1:00 9.45 1:30 10.28 2:00 10.89 2:30 11.37 3:00 11.76 3:30 12.09 4:00 12.36 4:30 12.59 5:00 12.79 5:30 12.97 6:00 13.13 6:30 13.27 7:00 13.40 7:30 13.52 8:00 13.64 8:30 13.75 9:00 13.86 9:30 13.96 10:00 14.06 10:30 14.16 11:00 14.26 11:30 14.36 12:00 14.46 （6） 温度变化图
5. 一厚度为 250mm 无限大平壁， 其导热系数λ=43+0.08t w/(m· k)， 平壁一侧温度为 250℃， 另一侧温度为 46℃， 试用数值方法确定平壁内的温度分布， 并确定通过该平壁的热流密度。 解： （1） 建立模型 本题属于非常物性，无热源的一维稳态导热过程 将平壁沿厚度方向（x 方向）划分为 N 个均匀相等的间距。 节点布置如图所示
输入 n, e, k, t1, tn+1, ti
N=10, t1=250, tn+1=46, ti=1, e=0.01, k=10000
迭代次数 k=0
tti=ti
No
k=k+1
Yes No k>m Yes
打印 m, ti, q
打印“不收敛”
停机
（4） 编写程序 本题使用 matlab 软件，所编写的程序如下： clear;clc; t=ones(1,11); % 设定各项初始值 q=ones(1,11); t(1)=250; t(11)=46; e=0.01; k=1; while 1 % 迭代程序 tt=t; for i=2:10 a=43+0.08*(t(i-1)+t(i))/2; b=43+0.08*(t(i)+t(i+1))/2; t(i)=(a*t(i-1)+b*t(i+1))/(a+b); end k=k+1; ttt=abs(t(5)-tt(5)); if(ttt<e) break; end end for i=2:11 % 计算热流密度 a=43+0.08*(t(i-1)+t(i))/2; q(i)=q(i-1)+a*(t(i)-t(i-1)); end k t q=q(i)/0.25 w=0:25:250;v=w/25+1;y=t(v); %绘制温度分布图 plot(w,y),xlabel('位置（mm）'),ylabel('温度（℃）') legend('平壁内的温度分布',0),grid （5） 计算结果 迭代次数 k=75 平壁温度分布见下表