计算传热学程序设计
中国石油大学(华东)
储建学院热能与动力工程系
《计算传热学程序设计》
设计报告
1引言
有关墙体传热量计算的方法是随着人们对房间负荷计算精度要求的不断提高而不断发展的.考虑辐射强度和周围空气温度综合作用，当外界温度发生周期性的变化时，屋顶内部的温度和热流密度也会发生周期性的变化。
计算题目
有一个用砖墙砌成的长方形截面的冷空气通道，其截面尺寸如图1所示。假设在垂直于纸面方向上冷空气及砖墙的温度变化相对较小，可近似地予以忽略。试计算稳态时砖墙截面的温度分布及垂直于纸面方向1米长度的冷量损失。设砖墙的导热系数为(m·℃)。内、外壁面均为第三类边界条件，外壁面：t
=30℃,h1=10W(m2·℃);内壁面：t f2=10℃, h2=4W(m2· ℃)。
f1
图1 砖墙截面
已知参数
砖墙的基本尺寸，砖墙的导热系数，外壁面的表面传热系数，对应的流体温度，内壁面的表面传热系数，对应的流体温度。
2 物理与数学模型
物理模型
由题知垂直于纸面方向上冷空气及砖墙的温度变化相对较小，可近似予以忽略，墙面为常物性，可以假设：
1）砖墙在垂直于纸面方向上没有导热。
2）由于系统是几何形状与边界条件是对称的，它的中心对称面就是一个绝热边界，这时只需求解1/4个对称区域就可以得到整个区域的解。 数学模型
考虑到对称性，取右下的1/4为研究对象，建立如图2的坐标系。
a
图2 砖墙的稳态导热计算区域
由上述的物理模型与上面的坐标系，该问题的数学模型可直接由导热微分方程简化而来，即
2222
0T 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)
1 1.1
1.1
(y y T h T
y
λ
==?-=-
? (4)
22(0.5,00.6)(0.5,00.6)
()f x y x y T h T T x λ
=<<=<
22(0.6,0.5 1.5)(0.6,0.5 1.5)
()f y x y x T h T T x
λ
=<<=<
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）