计算传热学程序设计

中国石油大学(华东)储建学院热能与动力工程系《计算传热学程序设计》设计报告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 Tx y ∂∂+=∂∂ (1)相应的边界条件是：1.10y T y =∂=∂1.50x T x=∂=∂ (2)110()f x x T h TT xλ==∂-=-∂ (3)111.11.1()f y y T h TT 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 在任意时刻均满足控制方程：[()()]0PP 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 sT 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 eya x δδλ=，()/W w wya x δδλ=，()/N n nxa y δδλ=，()/S s sxa 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.11y M δ=- (13)由题意知bc,cd,cf,af 均为第三类边界条件，ab,cg,de 为绝热边界。

可以采用元体能量平衡法得到各边界及拐点的差分方程（具体的边界节点及其拐点的离散过程见附录），分别为： 边界ab:22222W E NB T T Lr T T Lr++=+ (14) 拐点b ：22222W N fB T Lr T h yLrT T Lr h Lr yλλδλλδ++=++ (15)边界bc ：222222()2222W N S f B T Lr T T h LrT yT Lr h Lr yλλδλλδ+++=++ (16)拐点c ：22222222()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 ：22222()22222W E N fB 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 ：22222W N fB T Lr T h Lr xT T Lr h Lr xλλδλλδ++=++ (20) 边界de ：222()22W N S B T Lr T T T Lr++=+ (21)拐点e ：12121W S fB T Lr T h Lr xT T Lr h Lr xλλδλλδ++=++ (22) 上边界ef ：12121()22222E W S f B T T Lr T h LrT xT Lr h Lr xλλδλλδ+++=++ (23) 左边界af ： 121212()2222E S f B T Lr T TN h LrT yT Lr h Lr yλλδλλδ+++=++ (24)拐点a ：12121E N fB T Lr T h Lr yT T Lr h Lr yλλδλλδ++=++ (25) 拐点f ：12121()()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 可以得到以下结论：由图像可知所得到的温度曲线时连续的，并且在整个矩形界面的温度曲线应该是封闭的。

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

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

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

节点数的不同，可带来以下的区别：当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开始给相关参数赋值假设初始温度场 Gauss-Seidel 法求取温度场计算两次迭代误差是否满足计算精度输出至屏幕NY将温度分布输出至文件计算内外边界的热量并输出结束图5b图5c图5d图5e图5 计算区域的等温线图6结论有关墙体传热量计算的方法是随着人们对负荷计算精度要求的不断提高而不断发展的。

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

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

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

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

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

7参考文献【1】黄善波刘中良编著.计算传热学基础.中国石油大学（华东），2009杨世铭陶文铨编著.传热学（第四版）.高等教育出版社，2006刘衍聪编著.CAD技术基础.中国石油大学（华东），20068 附录附录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。