计算传热学程序设计

参考文献
[1] 黄善波 刘中良编着 计算传热学基础 [2] 杨世铭 陶文铨编着 传热学（第四版）北京高等教育出版社 [3] 王元明编 数学物理方程与特殊函数（第三版） 高等教育出版社
附录
1 数学模型的离散过程推导
2 p
x2
P
1 x
p x e
p
x
w
1 x
p E
x e
1 x
者吻合的比较好，可以说明所编制的数值解法的程序是正确的。相对误差先增大后减小，
增大的原因是此时温度接近零度，相对误差的基数比较小，所以造成相对误差较大，但
是此时的绝对误差并不大，在合理范围内，所以除去个别点外，都满足误差小于百分之
1。可以验证所编数值解法的程序是正确的。
（2） 空间步长对墙内壁的温度影响如图 6 及表 2。 在程序编写过程中用网格节点
图 12 墙壁厚度对内壁温度的影响
表 8 墙壁厚度对内壁温度的影响
时间（h）
0
15
15
15
15
10
6 结论
(1)整个墙壁经历了以下变化过程：首先外壁直接与室外冷空气接触，温度变化很 快，随着时间的推移，墙内各点的温度也开始变化，并影响到右边内墙壁的温度，慢慢 降低。
(2)对于墙内壁的温度随时间的变化，变化趋势总是由快到慢，最后重新达到稳态。 当改变墙壁厚度的时候，随着墙壁厚度的增加，对于外界温度的感应是越来越慢。这对 于一些对温度有要求的地反很重要，根据温度的变化作出适当的调整措施。
中国石油大学(华东)
储运与建筑工程学院热能与动力工程系
《计算传热学程序设计》 设计报告
学生姓名: 学 号： 专业班级： 指导教师
2012 年 7 月 7 日
1、设计题目
有一房屋的砖墙厚 δ= m ，λ= W/(m·℃)，ρc=×106 J/( m3·K)，室内温度 Tf1 保持 20℃不变，表面传热系数 h1=6W/(m2·℃)。开始时墙的温度处于稳定状态，内墙表 面温度 Tw1 为 15℃寒潮入侵后，室外温度 Tf2 下降为-10℃，外墙的表面传热系数为 35W/(m2·℃)。试分析寒潮入侵后多少时间内墙壁面方可感受到外界气温的变化。
x
a
0 p
2
T00
（18）
(
a
0 p
2
h1
)TN
x
1
TN
x
2 h1Tf 1
a0p 2
T0 N 1
（19）
2 程序清单
序号 程序名称
（1） 设计程序
程序功能
用数值方法求解 温度
计算结果 ① 表 1、4，图 5、8 ② 改变空间步长：表 2，图 6 ③ 改变时间步长：表 3，图 7 ④ 改变导热系数：表 5，图 9 ⑤ 改变墙外换热系数：表 6，图 10 ⑥ 改变墙内换热系数：表 7，图 11 ⑦ 改变墙壁厚度：表 8，图 12
程序特点：该程序有很强的适应性，一维常物性非稳态平壁导热问题都可以使用此 程序，只要适当更改边值条件即可。还可以进行修改解决非常物性问题。 程序中对输出节点，最大输出量都进行了控制，对计算结果的分析有很大帮助。而且 Thoms 算法的优点需要内存小，工作量小，程序设计简单。
程序流程图：首先对变量赋值，然后由初始条件建立初始温度场，接着从左边界， 内部节点，到右边界进行迭代，直到满足精度要求为止，最后输出结果，程序结束。程 序流程如下图 3。
分析解（℃）
15
数值解（℃）
15
相对误差（%）
0
5．00
图 5 平壁中心不同时刻的数值解和分析解
时间(h)
图 6 空间步长对墙内壁温度影响
表 2 空间步长对温度影响数据
N=51
N=101
N=201
（3） 时间步长对温度的影响如图 7 和表 3，根据图中曲线可以看出时间步长选择 50s，100s，200s 时基本重合，对墙内壁温度影响不大。
过建模，方程离散化，最终通过程序求解方程，得到图 8 和表 4。由图可以看出，开始
阶段，内墙壁温不变，随着时间的进一步深入，内壁温度开始降低，当很长时间后，温
度变化基本趋于平缓，直到再次平衡。根据图 8 就可以得到墙内壁温度开始发生变化的
时间。
时间（h)
表 4 墙内壁温度随时间变化数据表
温度（℃） 15 15
室外 寒流入侵
室内
0
x
图 1 墙壁简化图
已知参数 壁厚，墙壁导热系数，密度与比热容的乘积，室内和寒潮入侵后室外空气温度，室
内空气和外墙的表面传热系数，开始时稳定状态下的内墙表面温度。 求解
寒潮入侵多少时间后内墙壁面可感受到外界气温的变化？
2 物理与数学模型
物理模型 该墙面为常物性，可以假设：（1）其为无限大平面，（2）只有在厚度方向传热，
程序构成和方法:程序由主程序和一个子程序构成。主程序进行变量定义和各已知 参数的输入，以及左右边界节点和内部节点控制方程的输入；子程序 tdma 实现追赶法 用来计算每个节点新的温度。Thomas 算法求解过程分为两步:消元和回代。消元是从系 数矩阵的第二行起，逐一将每一行的非零元素消去一个，使原来的三元方程化为二元方 程。消元进行到最后一行时，二元方程就化为一元方程，直接得到最后一个未知数的值。 然后逐一往前回代，由各二元方程求出其它未知解。
Baidu Nhomakorabea
数对空间步长进行控制，为了观察空间步长对墙内壁温度的影响，表中选择了三个不同
的空间步长，分别为选取 51，101，201 个网格节点，则相应的空间步长为，， 。根据
不同步长时温度的变化曲线可以看出，空间步长对内墙壁的影响不大，当空间步长控制
在合理范围时可以忽略空间步长的影响。
时间（h）
0
表 1 分析解与数值解比较
左侧： h Tf
右侧： h=25 W/(m2·℃) Tf=-10℃
-
0
δx
δ
图 4 一维导热简化模型
直接根据公式得到解析解如下：
(,
0.
)
Cn
n1
exp(n2Fo) cos(n)
（10）
式中， Fo
a 2
,
x
，系数 Cn
应该使上述无穷级数在
0 是满足初始条件，由傅里叶
级数理论可得：
Cn
n
2 sin n cos n sin
e
1
x w
p P
Wp
x w
（1）
按照 Taylor 级数展开法，温度对时间的偏导有向前差分格式，中心差分格式和向
后差分格式，使用向后差分格式。
向后差分格式: 所以
由 x x 得
p
P
Pp
p1 P
Pp
p 1 P
a x
p E
x e
1
x e
1
x w
Pp
Wp
x w
开始
进行参数赋 值
用迭代法 求解温度 场
N
计数并判断是否完成 温度场计算
输入至屏幕 输入至文件 程序结束
引入过余温度：
Y
图 3 程序流程图
a c
T (x, 0) T0 (0 x )
T (x, ) x
x0
0
h
T
(
,
)
T
T
(x, x
)
x
T (x, ) T
（6） （7） （8）
（9）
T0=15℃
（3） （4）
（5） （6） （7）
如图 1： 其中
为常数
N-2
N-1 x
Δx
图 1 右边界的能量守恒
B N 2 Es
B
hr (Tfr
TN 1) ， N 2
N 2 N 1 x
N 2 N 1 x
hr (Tfr
TN1)
cx
N 1 0N 1
（8） （9） （10）
ac , x x 2
时间(h) 0
图 7 时间步长对温度的影响
表 3 时间步长对温度的影响
T=50（s） 15 15
时间(h) 0
T=100（s)
时间（h） 0
T=200（s） 15 15
2
4
6
8
图 8 墙内壁温度随时间的变化曲线
5 计算结果与分析
墙内壁温度分析
根据题目中要求，计算寒潮入侵多长时间后内墙壁可以感受到外界气温的变化，通
n
（11）
n 是超越方程的根，称为特征根。
其中 Bi h 。
程序验证
tan
n
Bi n
,n
1, 2, …
（12）
（1） 由模型可以得到相关信息然后进行编程，同等时间下计算出中心处温度的解
析解和数值解进行比较，数据记录在表 1。然后计算出相对误差，作图 5，观察数值解
与分析解的比较曲线。
由图表中可以发现，平壁中心不同时刻温度值的分析解和数值解相差不是很大，二
图 11 墙内表面传热系数对温度影响
表 7 墙内表面传热系数对温度影响
时刻
h=2
h=6
h=10
（h） W/(m2·℃) W/(m2·℃) W/(m2·℃)
墙壁厚度 δ 对内壁温度的影响 改变墙壁的厚度 δ，内壁温度将发生变化。在表 8 和图 12 中可以发现，不同厚度
的墙壁对外界温度的感应快慢是不一样的，随着墙壁厚度的增加，感应到外界温度变化 的时间越来越长，且温度变化越来越慢，墙壁越厚要越长的时间才能达到新的平衡
将式（12-b）(12-c)入式（12-a）以得到:
（11）
c x 2
x
hr
N
1
x
N 2
hrTfr
c x 2
T0 N 1
同理可以得到左边界点离散方程:
因此，离散方程为:
(hl
x
c x 2
)T0
hlTfl
x T1
c x 2
T00
（12） （13）
式中，
aPTP aETE aWTW b
导热系数，由于墙壁导热系数不变，故都等于 λ，△τ 为时间步长。由元体能量平衡法
可以得知左右边界节点的离散方程分别为：
左边界节点：
右边界节点：
（3）
（4）
离散方程的详细推导过程见附录。 程序设计
由物理模型可以知道本问题为一维导热问题，一维导热问题的离散方程在取遍所有 节点后形成的是三对角的代数方程组，采用追赶法进行求解。
此例中墙壁导热系数为常值，无源项。则可采用有限体积法对控制方程离散化，得 到离散方程为：
式中：
apTp aETE aWTW b
（2a）
aP
aE
aW
a
0 P
（2b）
aE
x
， aW
x
，
a
0 P
cx
（2c）
b a0pTp0
（2d）
其中的上标“0”表示此为上一时刻的值,分别为节点所在控制容积左右边界上的
表 6 墙外换热系数对温度影响
h=50
时刻
h=20
h=35
W/(m2·℃
（h） W/(m2·℃) W/(m2·℃) )
0
15
15
15
10
墙内壁传热热系数的影响 由图 11 可以看出，墙内表面传热系数对内表面温度影响较大，当传热系数比较小
时受墙外流体温度影响明显，传热系数越大受墙外流体温度影响越小，当到了极限大的 情况下，内表面温度则等于墙内流体温度，不再受墙外流体温度影响。
（14）
aP aE aW aP0
aE
x
，
aW
x
，
aP0
c x
b aP0TP0
式中，上标“0”表示上一时刻值， 为热扩散系数， 为时间步长。
由元体能量平衡法可以分别求出左右边界节点的离散方程，
左边界为:
（15） （16） （17）
右边界为:
(
a
0 p
2
h2
)T 0
x
T 1 h2Tf 2
时间（h)
温度（℃）
时间（h)
温度（℃）
导热系数 λ 对墙内壁温度的影响
墙的导热系数对内表面的影响，在图 9 和表 5 中发现，导热系数对内壁温度影响比
较大，λ=时，温度下降的趋势会更快，要比 λ=时下降快的多，下降速度更快，更短时
间内达到稳态，因此导热系数越大温度扩散越快，导热系数越小则温度变化越慢，需要
更长时间到达稳态，但是这时对于要求恒温的空间有好处，受波动影响更小。
表 5 导热系数对墙内壁温度的影响
时刻（h） 0
λ= 15 15
λ= 15 15
λ= 15 15
5
图 9 导热系数对墙内壁温度的影响
墙外换热系数 h 的影响 墙外表面传热系数对温度分布的影响，如图 10 和表 6 影响不大。
图 10 墙外换热系数对温度影响
左边界： 右边界：
3 数值处理与程序设计
T X
x0 h2 (Tx 0 Tf 2)
T X
x h1(Tx Tf 1)
（1c） （1d）
数值处理 采用外点法用均匀网格对求解区域进行离散化，得到的网格系统如图 2 所示。一共
使用了 0~N-1 共 N 个节点。 节点间距 δx 为:
图 2 墙壁内的网格划分
x
a
Pp
p1 P
Ep
x e
1
x e
1
x w
Pp
Wp
x w
x
2 x
P
x
E
x
W
c
x
0P
aPP aEE aW W b
aP aE aW aP0
aE
x
，
aW
x
，
aP0
c x
边界条件处理：右边界根据元体能量平衡法处理：
b aP0TP0
Φw
ΦB
（2）
没有纵向传热，则该问题转化为一维常物性无限大平面非稳态导热问题。 数学模型
以墙外表面为坐标原点，沿厚度方向为坐标正方向，建立坐标系。基于上述模型， 取其在 x 方向上的微元作为研究对象，则该问题的数学模型可描述如下：
c T ( T) x x
（1a）
初始条件：
（1b）
在两侧相应的边界条件是第三类边界条件，分别由傅立叶定律可描述如下：
4、模型与程序验证
模型
本题简化为厚度为 2 =的一维非稳态模型如图 4 所示，初始温度为 15℃，在其中 间建立坐标系，左两边为对流换热，且换热系数相同都为 h=25 W/(m2·℃)，且流体温 度 Tf=-10℃对于 x 0，列出其导热微分方程式及定解条件：
T
a
2T x2
(0
x
,
0)
（5）