计算传热学程序设计

中国石油大学(华东)
储建学院热能与动力工程系
《计算传热学程序设计》
设计报告
1引言
有关墙体传热量的方法是随着人们对房间负荷计算精度要求的不断提高而不断的.考虑辐射强度和周围空气温度综合作用，当外界温度发生周期性的变化时，屋顶内部的温度和热流密度也会发生周期性的变化。

计算题目
有一个用砖墙砌成的长方形截面的冷空气通道，其截面尺寸如图1所示。

假设在垂直于纸面方向上冷空气及砖墙的温度变化相对较小，可近似地予以忽略。

试计算稳态时砖墙截面的温度分布及垂直于纸面方向1米长度的冷量损失。

设砖墙的导热系数为(m·℃)。

内、外壁面均为第三类边界条件，外壁面：t f1=30℃,h1=10W(m2·℃);内壁面：t f2=10℃, h2=4W(m2·℃)。

图1 砖墙截面
已知参数
砖墙的基本尺寸，砖墙的导热系数，外壁面的表面传热系数，对应的流体温度，内壁面的表面传热系数，对应的流体温度。

2 物理与数学模型
物理模型
由题知垂直于纸面方向上冷空气及砖墙的温度变化相对较小，可近似予以忽略，墙面为常物性，可以假设：
1）砖墙在垂直于纸面方向上没有导热。

2）由于系统是几何形状与边界条件是对称的，它的中心对称面就是一个绝热边界，这时只需求解1/4个对称区域就可以得到整个区域的解。

数学模型
考虑到对称性，取右下的1/4为研究对象，建立如图2的坐标系。

a
图2 砖墙的稳态导热计算区域
由上述的物理模型与上面的坐标系，该问题的数学模型可直接由导热微分方程简化而来，即
22220T T x y ∂∂+=∂∂ (1)
相应的边界条件是：
1.1
0y T y =∂=∂
1.5
0x T x
=∂=∂ (2)
110
()f x x T h T
T x
λ
==∂-=-∂ (3)
111.1
1.1
()f y y T h T
T y
λ
==∂-=-∂ (4)
22(0.5,00.6)(0.5,00.6)
()f x y x y T h T T x λ
=<<=<<∂-=-∂ (5)
22(0.6,0.5 1.5)(0.6,0.5 1.5)
()f y x y x T h T T x
λ
=<<=<<∂-=-∂ (6)
3数学模型的离散化
采用外点法对求解区域进行离散化，其中ab 方向上取N 1个节点，af 边界上取M 个节点，bc 边界取M 1个节点，cd 边界取N 2个节点，de 边界取M 2个节点，ef 边界取N 个节点，离散后的求解区域如图3所示。

图3 求解区域的离散化
采用Taylor 级数展开法得到内部节点的差分方程，求解区域内任一节点P 在任意时刻均满足控制方程：
[
()()]0P
P T T x x y y
λλ∂∂∂∂+=∂∂∂∂ (7) Taylor 级数展开法的思想是用差商表达式代替方程中的各阶导数。

按照一维问题的处理方法，将方程中的各阶导数的差商表达式代入到上式，有
11
[
()][()()][]()()P W P E P P e w e w e w
T T T T T T T x x x x x x x x λλλλλδδδδ--∂∂∂∂=-=-∂∂∂∂
(8)
11
[
()][()()][]()()P N P P S P n s n s n s
T T T T T T T y y x y y y x x λλλλλδδδδ--∂∂∂∂=-=-∂∂∂∂ (9) 将以上(8)，(9)代入 (7)式可得到
P P E E W W N N S S a T a T a T a T a T =+++ (10)
式中
()/E e e
y
a x δδλ=
，()/W w w
y
a x δδλ=
，()/N n n
x
a y δδλ=
，()/S s s
x
a y δδλ=
(11)
可得：
222()/[2(1)]p W E r S r N r T T T L T L T L =++++ (12)
式中，
r x L y δδ=
， 1.51x N δ=-， 1.1
1
y M δ=- (13)
由题意知bc,cd,cf,af 均为第三类边界条件，ab,cg,de 为绝热边界。

可以采用元体能量平衡法得到各边界及拐点的差分方程（具体的边界节点及其拐点的离散过程见附录），分别为： 边界ab:
22
222W E N
B T T Lr T T Lr ++=+ (14)
拐点b ：
2
2222W N f
B T Lr T h yLrT T Lr h Lr y
λλδλλδ++=++ (15)
边界bc ：
222222()2222W N S f B T Lr T T h LrT y
T Lr h Lr y
λλδλλδ+++=
++ (16)
拐点c ：
2
2222222()33()W N S E f B T Lr T Lr T T h x y LrT T Lr h x y Lr
λλλλδδλλδδ+++++=
+++ (17)
边界cd ：
2
2222()22222W E N f
B T T Lr T h xLrT T Lr h Lr x
λλδλλδ+++=++ (18)
虚拟边界cg ：
22
()22W E N S B T T Lr T T T Lr
+++=+ (19) 拐点d ：
2
2222W N f
B T Lr T h Lr xT T Lr h Lr x
λλδλλδ++=
++ (20) 边界de ：
22
2()
22W N S B T Lr T T T Lr
++=+ (21) 拐点e ：
1
2121W S f
B T Lr T h Lr xT T Lr h Lr x
λλδλλδ++=
++ (22) 上边界ef ：
1
2121()22222E W S f B T T Lr T h LrT x T Lr h Lr x
λλδλλδ+++=
++ (23) 左边界af ：
1212
12()2222E S f B T Lr T TN h LrT y
T Lr h Lr y
λλδλλδ+++=
++ (24)
拐点a ：
1
2121E N f
B T Lr T h Lr yT T Lr h Lr y
λλδλλδ++=
++ (25) 拐点f ：
1
2121()()
E S f B T Lr T h Lr x y T T Lr h Lr y x λλδδλλδδ+++=+++ (26)
4程序编写及验证
程序设计的思路
首先 ，对计算区域进行均匀格划分，给出第一部分的x 方向的节点数，计算出第一部分的y 轴方向的节点数以及第二部分x ，y 方向的节点数，并计算出整个计算区域的x ，y 方向的总节点数。

然后对其温度场的假设，在开始编程时将温度场划分为两部分，但是在运用Tecplot 软件对计算区域进行绘图时不能对计算区域的温度场划分，所以在假设时将其划为一个区域。

再下一个环节是由Gauss-Seidel 迭代法计算各结点温度，结点编程即为计算区域边界点和拐点的编程，顺序是y 轴方向上由下到上，x 方向上由左到右，有关内部节点的编程。

接下来是对计算两次迭代间的最大误差，判断是否满足计算精度，输出计算结果，等温线的数据文件的输出的编程。

最后就是对编程的校核，即为冷量的计算输出的编程。

程序流程图（如图4）
程序的验证
1) 验证的必要性说明
计算传热学的验证是经过科学设计的程序，或当问题有精确地理论解时，用验证程序去说明程序的准确性与可行性，或对照理论解检验你所处理的问题的正确性。

2)本例的验证过程
程序的验证是通过检验内外边界所传出的热量相等来验证。

在本题中所研究的对象在几何形状和边界条件是对称的，所以取了1/4的单元来研究。

当所求的温度分布正确时，可以得到内外边界的热量是相等的。

通过编程可以知道外边界传出的热量Q1,Q2，存在的误差在可以接受的范围之内，即证明了程序的可靠性。

5 计算结果和分析
通过程序结果可以知道外边界传出的热量Q1，Q2（如下一页图5）。

可知Q1基本等于Q2，存在的误差。

误差存在的原因可能由于计算过程和编程过程中精确度的取舍有关。

由图5 可以得到以下结论：
(1)由图像可知所得到的温度曲线时连续的，并且在整个矩形界面的温度曲
线应该是封闭的。

(2)截面上的温度曲线是光滑的，没有特别的凸起和凹陷，说明其温度具有
一定规律的分布。

(3)温度曲线没有交叉，且在砖墙的内壁面附近温度曲线密集，在外壁面附
近温度曲线较稀疏。

通过不同的取节点数，可以得到不同的冷量，从而得到不同的误差。

节点数的不同，可带来以下的区别：
当N=16，M=12时，Q1=，Q2=，ε1=。

(如图5a)
当N=31，M=23时，Q1=，Q2=，ε2=。

(如图5b)
当N=61，M=45时，Q1=，Q2=，ε3=。

(如图5c)
当N=166，M=122时，Q1=，Q2=，ε4=。

(如图5d)
由上数据和图5a-图5-d可知，所取节点数越多，靠近内壁附近的温度曲线越密集，墙壁外壁附近的温度曲线分布越疏松，但是，内部导热系数较小，外部导热系数较大，得到结果的误差越大。

改变墙体的厚度，可以得到不同的冷量，从而得到不同的误差。

厚度的不同可以带来以下误差：
当厚度为，N=17，M=13时，Q1=，Q2=，ε5=0。

(如图5e)
当厚度为，N=16，M=12时，Q1=，Q2=，ε1=。

(如图5a)
当厚度为，N=15，M=11，Q1=，Q2=，ε6=0。

由以上数据和图5a，图5e可知，墙体厚度对墙体内温度曲线的分布影响很小，冷量改变基本不变。

图4 程序流程图
图5a
图5b
图5c
图5d
图5e
图5 计算区域的等温线图
6结论
有关墙体传热量的方法是随着人们对负荷计算精度要求的不断提高而不断的。

本次设计的题目中墙体的截面点平面结构为矩形，可能存在拐点处的传热的不平衡，从而导致加热或者冷却的传热效率差的现象。

如果将其平面结构改成圆形，或许会更优越一点。

有上述结果分析可得到以下结论：
(1)取节点数越多误差越大。

(2)墙体的厚度对传热误差不大。

(3) 通过对温度曲线的观察可以看出温度的升高或降低并不是线性的，越靠近内壁面温降越大，靠近外壁面的温降比较平滑，温降也比较小。

7参考文献
【1】黄善波刘中良编著.计算传热学基础.中国石油大学（华东），2009【2】杨世铭陶文铨编著.传热学（第四版）.高等教育出版社，2006
【3】刘衍聪编著.CAD技术基础.中国石油大学（华东），2006
8 附录
附录1 附件中程序清单
1.许生举-01………主程序(N=16，M=12厚度为如图5a。

2. 许生举-02………修改网格数目后的程序(N=31，M=23厚度为如图5b。

3. 许生举-03………修改网格数目后的程序(N=61，M=45厚度为如图5c。

4. 许生举-04………修改网格数目后的程序(N=166，M=122厚度为如图5d。

5. 许生举-05………修改墙体的厚度( N=17，M=13 壁厚为如图5e。

附录2 边界的温度推导
各个边界节点及拐点温度的推导
在第三类边界条件中，由于边界上的温度未知，因此为了构成封闭的方程组，必须补充边界节点的离散方程。

这里采用了元体能量平衡法导出各边界的离散方程。

1）第一部分下边界ab
图 6 ab的元体能量平衡法
对如图6所示的边界节点，考虑能量平衡。

根据能量守恒，有
Φ+Φ+Φ+Φ+Φ=(27)
W E N B V
无内热源，不用考虑内热源Φv的影响。

ΦB为由边界流入控制容积的热量，规定流入边界的为正，但是ab边界为绝热边界即
Φ=(28)
B
式(27)中的ΦW，ΦN，ΦE为通过控制容积界面导入微元体的热量。

假设节点间的温度分布按线性分布，根据Fourier定律有：
y
2()W B
W W
T T x λδ-∆Φ= (29)
2()E B
E E
T T y x λδ-∆Φ=
(30)
()N B
N N
T T
x y λδ-Φ=∆
(31)
将(27)，(28)，(28)，(30)，(31)联立，可解得下边界ab 节点的差分方程为：
22
222W E N
B T T Lr T T Lr ++=+ (32 )
2) 第一部分右下拐点b
图7 b 的元体能量平衡法
对如图7所示的边界节点，考虑能量平衡。

根据能量守恒，有
0W R N r V Φ+Φ+Φ+Φ+Φ= (33)
无内热源，不用考虑内热源Φv 的影响，同时ab 为绝热边界，则
0r Φ=
(34)
ΦB 为由边界流入控制容积的热量，规定流入边界的为正，即
222()()2
B f B x h T T ∆Φ=- (35) 式(33)中的ΦW ，ΦN 为通过控制容积界面导入微元体的热量。

假设节点间的温度分布按线性分布，根据Fourier 定律有：
2()W B
W W T T y x λδ-∆Φ=
(36)
2()N B
N N T T x y λδ-∆Φ=
（37）
将(33)，(34)，(35)，(36)，(37)联立，可解得第一部分右下拐点b 的差分方程为：
2
2222W N f B T Lr T h yLrT T Lr h Lr y λλδλλδ++=++ (38)
图8 bc 的元体能量平衡法
3) 第一部分右边界bc
对如图8所示的边界节点，考虑能量平衡，根据能量守恒，有
0W S N B V Φ+Φ+Φ+Φ+Φ= (39)
无内热源，不用考虑内热源Φv 的影响， ΦB 为由边界流入控制容积的热量，规定流入边界的为正，即
222()B f B h y T T Φ=∆- (40)
式(39)中的ΦW ，ΦS ，ΦN 为通过控制容积界面导入微元体的热量。

假设节点间的温度分布按线性分布，根据Fourier 定律有：
()W B W W
T T y x λδ-Φ=∆ (41) 2()S B S S T T x y λ
δ-∆Φ= (42)
2()N B N N
T T x y λ
δ-∆Φ= (43) 将(39)，(40)，(41)，(42)，(43)联立，可解得第一部分右边界bc 的差分方程为：
222222()2222W N S f B T Lr T T h LrT y
T Lr h Lr y λλδλλδ+++=++ (44)
图9 c 的元体能量平衡法
4) 第一部分右上拐点c
对如图9所示的拐点，考虑能量平衡。

根据能量守恒，有
0W B N S V E b Φ+Φ+Φ+Φ+Φ+Φ+Φ= (45)
无内热源，不用考虑内热源Φv 的影响，Φb 为由边界流入控制容积的热量，规定流入边界的为正，即
222()()2
b f B y h T T ∆Φ=- (46) ΦB 为由边界流入控制容积的热量，规定流入边界的为正，即
222()()2
B f B x h T T ∆Φ=- (47) 式(45)中的ΦW ，ΦS ，ΦN ，ΦE 为通过控制容积界面导入微元体的热量。

假设节点间的温度分布按线性分布，根据Fourier 定律有：
()W B W W
T T y x λδ-Φ=∆ (48) 2()2()S B S S
T T x y λδ-∆Φ= (49) ()N B N N T T x
y λδ-Φ=∆ (50) 2()2()E B E E
T T y y λδ-∆Φ= (51) 将(45)，(46)，(47)，(48)，(49)，(47)，(51)联立，可解得拐点c 的差分方程为：
2
2222222()33()W N S E f B T Lr T Lr T T h x y LrT T Lr h x y Lr λλλλδδλλδδ+++++=+++
(
5
2
)
5) cd
对如图10所示的边界节点，考虑能量平衡，根据能量守恒，有
0W E N B V Φ+Φ+Φ+Φ+Φ= (53)
无内热源，不用考虑内热源Φv 的影响。

ΦB 为由边界流入控制容积的热量，规定流入边界的为正，即
222()B f B h x T T Φ=∆- (54)
式 (53) 中的ΦW ，ΦE ，ΦN 为通过控制容积界面导入微元体的热量。

假设节点间的温度分布按线性分布，根据Fourier 定律有：
2()W B W W
T T y x λδ-∆Φ= (55) 2()E B E E
T T y x λδ-∆Φ= (56)
()N B N N
T T x y λδ-Φ=∆ (57) 将(53)，(54)，(55)，(56)，(57)联立，可解得下边界cd 节点的差分方程为：
2
2222()22222W E N f B T T Lr T h xLrT T Lr h Lr x λλδλλδ+++=++ (58)
图10 cd 的元体能量平衡法
6) 第二部分右下拐点d
如图11所示的边界节点，考虑能量平衡，根据能量守恒，有
0W B N r V Φ+Φ+Φ+Φ+Φ= (59)
无内热源，不用考虑内热源Φv 的影响，同时de 为绝热边界，则
0r Φ=
(60)
ΦB 为由边界流入控制容积的热量，规定流入边界的为正，即
222()()2
B f B x h T T ∆Φ=- (61) 式 (59) 中的ΦW ，ΦN 为通过控制容积界面导入微元体的热量。

假设节点间的温度分布按线性分布，根据Fourier 定律有：
2()W B W W
T T y x λδ-∆Φ= (62) 2()N B N N
T T x y λδ-∆Φ= (63)
将(59)，(60)，(61)，(62)，(63)联立，可解得第二部分右下拐点d 的差分方程为:
2
2222W N f B T Lr T h Lr xT T Lr h Lr x λλδλλδ++=++ (64)
图11 d 的元体能量平衡法
7) 第二部分右边界 de
对如图12所示的边界节点，考虑能量平衡。

根据能量守恒，有
0W S N r V Φ+Φ+Φ+Φ+Φ=
(65)
无内热源，不用考虑内热源Φv 的影响, Φr 为由边界流入控制容积的热量，但是de 边界为绝热界面，即
0r Φ= (66)
式(65)中的ΦW ，ΦS ，ΦN 为通过控制容积界面导入微元体的热量。

假设节点间的温度分布按线性分布，根据Fourier 定律有：
()W B W W
T T y x λδ-Φ=∆ (67)
2()S B S S
T T x y λδ-∆Φ= (68) 2()N B N N T T x y λ
δ-∆Φ= (69) 将(66)，(67)，(68)，(69)，(65)联立，可解得第二部分右边界de 的差分方程为：
22
2()22W N S B T Lr T T T Lr ++=+ (70)
图12 de 的元体能量平衡法
8) 第二部分右上拐点e
根据计算图示可知离散方法与第二部分右下拐点相似，在此省略推导过程，则第二部分右上拐点e 的差分方程为：
1
2121W S f B T Lr T h Lr xT T Lr h Lr x λλδλλδ++=++ (71) 9) 第二部分上边界ef
对如图13的边界节点，考虑能量平衡，根据能量守恒，有
0W S E B V Φ+Φ+Φ+Φ+Φ= (72)
无内热源，不用考虑内热源Φv 的影响，ΦB 为由边界流入控制容积的热量，规定流入边界的为正，即
121()()B f B h x T T Φ=∆- (73)
式(72)中ΦE ，ΦS 为通过控制容积界面导入微元体的热量。

假设节点间的温度分布按线性分布，根据Fourier 定律有：
2()W B
W W
T T y x λ
δ-∆Φ= (74)
()S B
S S
T T x
y λδ-Φ=∆ (75)
图13 ef 的元体能量平衡法 2()E B
E E
T T y x λ
δ-∆Φ= (76)
将(72)，(73)，(74)，(75)，(76) 联立，可解得第二部分上边界ef 分方程为：
1
2121()22222E W S f B T T Lr T h LrT x T Lr h Lr x
λλδλλδ+++=
++ (77) 10) 左边界 af
对如图14的边界节点，考虑能量平衡，根据能量守恒，有
0B S N E V Φ+Φ+Φ+Φ+Φ= (78)
无内热源，不用考虑内热源Φv 的影响，ΦB 为由边界流入控制容积的热量，规定流入边界的为正，即
121()B f B h y T T Φ=∆- (79)
式(78)中 ΦE ，ΦS ，ΦN 为通过控制容积界面导入微元体的热量。

假设节点间的温
度分布按线性分布，根据Fourier 定律有：
()E B
E E
T T y
x λδ-Φ=∆ (80) 2()S B
S S
T T x y λ
δ-∆Φ= (81)
2()N B
N N
T T x y λ
δ-∆Φ= (82)
图14 af 的元体能量平衡法
将(78)，(79)，(80)，(81)，(82) 联立，可解得左边界的差分方程为：
1212
12()2222E S f B T Lr T TN h LrT y
T Lr h Lr y
λλδλλδ+++=
++ (83)
11) 第一部分左下拐点 a
拐点a 的性质与拐点b 的性质相似，所以离散过程在此省略。

则拐点a 的差分方程为：
1
2121E N f
B T Lr T h Lr yT T Lr h Lr y
λλδλλδ++=
++ (84) 12 ) 左上拐点f
对如图15所示的拐点，考虑能量平衡。

根据能量守恒，有
0E B S V b Φ+Φ+Φ+Φ+Φ= (85)
无内热源，不用考虑内热源Φv 的影响, Φb 为由边界流入控制容积的热量，规定流入边界的为正，即
12
1(
)()2
b f B x h T T ∆Φ=- (86)
ΦB 为由边界流入控制容积的热量，规定流入边界的为正，即
12
1(
)()2
B f B y h T T ∆Φ=- (87) 式(85)中ΦS ，ΦE 为通过控制容积界面导入微元体的热量。

假设节点间的温度分布按线性分布，根据Fourier 定律有：
2()S B
S S
T T x y λ
δ-∆Φ=
(88)
2()E B
E E
T T y y λ
δ-∆Φ=
(89)
将(85)，(86)，(87)，(88)，(89)联立，可解得拐点f 差分方程为：
1
2121()()
E S f
B T Lr T h Lr x y T T Lr h Lr y x λλδδλλδδ+++=+++ (90)
图15 f 的元体能量平衡法 13) 虚拟边界cg 和内部节点
对如图16所示的边界节点，考虑能量平衡。

根据能量守恒，有
W E N S V Φ+Φ+Φ+Φ+Φ= (91)
式(91)中ΦW ，ΦN ，ΦS ，ΦE 为通过控制容积界面导入微元体的热量。

假设节点间的温度分布按线性分布，根据Fourier 定律有：
()W B
W W
T T y
x λδ-Φ=∆ (92) ()N B
N N
T T x
y λδ-Φ=∆ (93)
()S B
S S
T T x
y λδ-Φ=∆ (94) ()E B
E E
T T y
x λδ-Φ=∆ (95) 将(91)，(92)，(93)，(94)，(95)联立，可解得内部节点及假想边界cg 的差分方程为：
22
()
22W E N S B T T Lr T T T Lr
+++=+ (96)
图16 内部节点及虚拟边界cg的元体能量平衡法
热能与动力工程系
《计算传热学程序设计》成绩考核表
教师签字：。